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Stable Artificial Dissipation Operators for 
Finite Volume Schemes on Unstructured Grids 

Magnus Svard* Jing Gong*and Jan Nordstrom* 


Abstract 

Our objective is to derive stable first-, second- and fourth-order ar- 
tificial dissipation operators for node based finite volume schemes. Of 
particular interest are general unstructured grids where the strength 
of the finite volume method is fully utilised. 

A commonly used finite volume approximation of the Laplacian 
will be the basis in the construction of the artificial dissipation. Both 
a homogeneous dissipation acting in all directions with equal strength 
and a modification that allows different amount of dissipation in dif- 
ferent directions are derived. Stability and accuracy of the new opera- 
tors are proved and the theoretical results are supported by numerical 
computations. 


1 Introduction 

In computational fluid dynamics, edge based finite volume (FV) approxima- 
tions are widely used (see [1, 2, 3, 4, 5, 6, 7, 8]). The main advantage of those 
schemes is a property called grid transparency by Haselbacher et al in [1]. 
Grid transparancy means that the algorithm only needs information about 
what nodes connect to each other and thus works equally well on structured 
as well as unstructured grids. This property is essential for efficiency and 
applicability of the scheme. 

* Center for Turbulence Research, Building 500, Stanford University, Stanford, CA 
94305-3035, USA, e-mail: svard@stanford.edu 

department of Information Technology, Uppsala University , Uppsala, Sweden, 
department of Computational Physics, Division of Systems Technology, The Swedish 
Defense Research Agency, SE-164 90 Stockholm, Sweden and Department of Information 
Technology, Uppsala University , Uppsala, Sweden. 
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Numerical computations often require artificial dissipation to remove un- 
physical high-frequency oscillations. Usually second- or fourth-order artificial 
dissipation are used. (The use of order refers to the order of accuracy.) When 
shocks are present first-order dissipation need to be used. (See for example 
[9] and [10].) The following properties are desirable. The artificial dissipation 
should: 

1. Reduce spurious oscillations. 

2. Preserve the order of accuracy of the numerical scheme. 

3. Preserve stability of the numerical scheme. 

Property 1 is achieved by using undivided differences. To preserve the accu- 
racy (Property 2) of the numerical scheme a sufficiently high-order derivative 
corresponding to the undivided difference need to be used. For a fourth-order 
accurate numerical scheme, a fourth-order undivided difference is added. One 
could also add a second-derivative approximation scaled with the grid size 
to obtain fourth-order accuracy. However, that implies that the damping of 
the highest frequency goes to zero as the grid is refined. With a fourth-order 
undivided difference the damping of the highest frequency is independent 
of the grid size. The treatment of Properties 1 and 2 is well-known and a 
variety of different dissipations have been proposed (See [9]). However, for 
unstructured finite volume schemes Property 3 have recieved little attention 
until now. We will focus on stability properties to obtain different artificial 
dissipation operators that satisfy all three properties even on unstructured 
grids. 

In [6], Nordstrom et al considered a first derivative approximation such 
that convective terms can be implemented in a stable and accurate man- 
ner. The stability proofs include boundary conditions since the scheme is 
interpreted in a summation-by-parts framework. This work was continued 
in [11] where an approximation of the Laplacian was interpreted in the same 
framework such that schemes consisting of first derivatives and Laplacians 
can be proven stable. In this paper we aim to construct first-, second- and 
fourth-order dissipation and in order not to destroy the stability of the main 
scheme, the artificial dissipation has to be bounded in the same norm as the 
main scheme. In [6] and [11] the norm is weighted with the finite volumes 
that discretise the domain. 

The contents of this report are divided as follows; in section 2 the general 
finite volume approximation is derived; in section 3 a second-order dissipa- 
tion is derived; section 4 contains a derivation of a fourth-order dissipation 
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operator; in section 5 the dissipation operator is modified to act with differ- 
ent strength in different directions. Section 6 shows numerical computations 
using the new dissipation operators and in section 7 conclusions are drawn. 

2 The Approximation of the Laplacian 

We aim to construct dissipation operators based on the application of the 
Laplacian. Therefore, we begin by deriving the standard node centred finite 
volume approximation of the Laplacian (see [1, 2, 3, 4, 5, 11] ). Since our main 
interest is to prove stability for the time dependent problem, we consider the 
heat equation, 

u t = A u. (1) 

Integration of (1) over the domain V] yields, 



where Gauss’ theorem is used. N denotes the outward pointing unit normal 
vector such that = u N = Vw • N. Further, let V t be an n-sided polygon 
with sides dsi n . 

Given any grid, let r t denote a grid point. With a slight abuse of notation 
we let Vi be defined as the measure of the volume inside the so called dual 
grid around r,. The dual grid is in turn defined as the straight lines drawn 
between the centres of mass of the cells with r, as a vertex and the midpoints 
of the edges from r,, see Figure 1. 



Figure 1: A generic 2D grid. Solid lines are the grid lines and dashed lines 
correspond to the dual grid. 
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Further, dsi n is defined as the sum of the length of the “centre of mass- 
midpoint-centre of mass” lines passing over one edge (see Figure 1). Let 
t hi = r, — r n \. Finally, let N t denote the set of indices of points being 
neighbours to r^. 

For an interior point, r t , an approximation of (2) is, 

Vi • (■ v t )i = ^2 — — ds in . (3) 

neNi Tm 


and for a boundary point b, 

Vb ■ ( Vt)b = ^2 — — ~ ds t>n + ( V N ) b ds bb , (4) 

^ AT F nb 

neN b 

where ds bb is the length between two midpoints of edges at the boundary. 
The scheme (3) and (4) can be summarised in matrix form as, 

Pv t = Q a v = {A + DS)v (5) 

where P is a diagonal matrix with Vi on the diagonal. DSv holds the terms 
(v N )bdsbb- Av represents the remaining terms (essentially the scheme for the 
interior points). v t and v are vectors with components (v t ) t and Vi respec- 
tively. In [11], A is proven to be a symmetric negative-definite matrix. 


3 Second-Order Artificial Dissipation 

A second-order dissipation operator is obtained by, 

u t = e 2 A u, 0 < x < 1, (6) 

where e is a small positive number (compare with the scaling h 2 in the discrete 
artificial dissipation). If we apply the energy method to (6) we obtain, 

/ uu t dx = e 2 uu x \l — e 2 / u 2 dx. (7) 

Jo Jo 

We see that eu = 0 or eu x = 0 will result in a well posed problem. (If e — > 0 
the boundary conditions vanish. This is the analogue of numerical boundary 
conditions.) In the discrete setting we have, 

v t = ^P^Qav. (8) 
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To analyse the discretisation of (6) we introduce the norm ||u|| 2 = v T Pv. 
Apply the energy method, 

|M| \ = v T Pv t + vJPv = h 2 2v T (A + DS)v = h 2 (2v T Av + 2v T (DSv)). (9) 

Since A is negative definite the discretisation is stable if v = 0 or DSv = 0 
at the boundary. 0 denotes a vector with all entries zero. With the SAT 
technique (Simultaneous Approximation Term, see [12, 13, 14]) we impose 
DSv = 0 weakly by adding, 

h 2 P~ x {DSv- 0 ), (10) 

to the right-hand side of (8). Finally, we have the artificial dissipation, 

v t = h 2 P~ x Av. (11) 

It is easy to determine the following sizes from (4). dim denotes the number 
of space dimensions. 


P-'QaV ~ 0(1), 

DSv ~ ds bb ~ 0(h dim '), 

p-‘~0( (12) 

dsn 

Av ~ Oik*™- 1 ). 

This leads to second-order interior accuracy and first-order boundary accu- 
racy, which result in globally second-order accuracy. 

Remark A scalar equation has been considered when deriving the artificial 
dissipation. That equation represents the treatment of a single variable in 
a system of partial differential equations such as the Euler or Navier-Stokes 
equations. The numerical boundary conditions discussed above, only describe 
how to close the scheme. They do not affect the physical boundary conditions 
and need not be changed if another equation is considered. Note that no 
boundary conditions are imposed in Equation (11) and that no boundary 
conditions should be imposed since it is an approximation of u t = 0. 

3.1 Some Remarks on First Order Dissipation 

It is well-kown that a shock capturing dissipation needs to be first order 
which can be obtained by dividing the second-order dissipation by the grid 
size. We propose one such scaling in Section 6 on general unstructured grids. 
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Another issue is the scaling of the dissipation for linear problems. To 
address this question we consider two different kinds of grid independence. 
First, if the entire problem is rescaled such that the grid size is doubled but 
with the same number of points (that is a twice as big grid) and the equations 
are rescaled accordingly so that the exact solution at each grid point is the 
same on the small and big grid, then that should be the case for the numerical 
solution as well. We illustrate this with an example. Consider, the periodic 
problem u t + au x = 0 with initial data u(x, 0) = f(x) on 0 < x < 1. Then 
the solution is u(x, t) = /(£) where £ = x — at+an and n is an integer chosen 
such that £ G [0, 1]. 

Next, consider the periodic problem v t +2 av y = 0 with v(y , 0) = f(y/ 2) = 
f(x) on 0 < y < 2. The solution is v(y , t ) = f(rj/2) where r] = y — 2 at + 2 an. 
n is an integer chosen such that rj G [0, 2], 

The two examples have identical solutions at the points located at the 
same ratio of the total interval. For example u(l/8,t) = v(l/4, t). This 
is easily understood since the two problems connects through the mapping 
y = 2x such that u x = u y y x = 2 u y . 

Next, we turn to the discrete problems on the same domains. Let u t + 
aD x u = A x u where D x is the discrete ^-derivative on x G [0, 1] with grid 
spacing h and A x is an artificial dissipation operator. Further, ii(x, t . 0) = 
f(xi). Let v t + 2 aD y v = A y v be the same discretisation with grid spacing 
2h on y G [0,2] and u(yi, 0) = f(xi). Note that we have the same number 
of grid points in the two problems. The two problems are indentical under 
the transformation y = 2x and so are the discrete solutions. Hence, at grid 
point i we have, u(xi,t) = v(yi,t)- Note that we have not assumed anything 
regarding the size of the dissipation so this similarity applies to both first- 
and second-order dissipation. 

The second property of importance is the scaling of the dissipation for a 
given linear problem and for different grid sizes. 

We consider the periodic problem u t + au x = 0on0<x<l described 
previously in this section. It is discretised by v t + aD x v = A x v where D x v 
is the standard second-order (non-dissipative) approximation. A x = h J A x 
where A x v is the undivided second-derivative approximation times a constant 
c. j = 1,2,3... yields a jth-order dissipation. Written in component form 
the N + 1-point discretisation is, 


/ \ . W+i W— l .j W+i 2 v n + v n —i 

( v„)t + a ^ = h J c , n = 1..N, v 0 = v N . 


2h 


h 2 
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Expand the solution in a Fourier series with h = 1/N, 


N - 1 

v(t) n = J2Ht)me ir ™ n i 

m— 0 


u n = 2irx n e [0, 2it(N - l)/N], 


N - 1 

v{t) m = ^2 v{x n )e~ %moJn h. 

n — 0 


( 13 ) 


Insert the Fourier expansion into the error equation, 

N - 1 

^(0(t) m )te® ,T1Wn — v{t) m e im,JJn {— ^zsin(m27rh) + 2c/i J_2 (cos(ra27r/i) — l))j = 0. 

m=0 

Note that with this choice of transform m = 0 corresponds to a constant and 
m = N — 1 is the least oscillating mode while m = N/2 is the most oscillating 
(usually called the 7r-mode since with m = N/2 we have m2-Kh = 7r). 

The purpose of the dissipation is to damp the unresolved modes while 
leaving the resolved modes undisturbed. Hence a low frequency mode should 
have less and less damping as the grid is refined and the highest mode should 
experience the same damping. Therefore, we study the behaviour of the mode 
m = N/2. The equation for that mode is, 

(v(t) N / 2 ) t = ^ (^) at /2 ( — ~?sin(27r/iiV/2) + 2ch- ?_2 (cos(27rhiV/2) — 1) = 

f V 

-v{t) N /2Ach 3 ~ 2 . (14) 

Now, it is easily seen that the only h independent solution is achieved with 
j = 2. If j = 1 the highest frequency will have increased damping as h 0 
and if j > 2 the damping will diminish with smaller h. 

These conclusions hold for any constant ratio N/x where x is a valid 
choice that gives an existing mode. Hence the choice j = 2 imposes the same 
damping, independent of h , on the portion of the modes that are unresolved. 
However, if a fixed mode N = Nf ix is considered the damping on that mode 
will decrease with h since it will be more and more resolved. 

The general conclusion for a pth-order scheme is that the undivided dif- 
ference corresponding to the pth derivative gives the grid independent dissi- 
pation. 
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4 Fourth- Order Artificial Dissipation 

The most obvious fourth-order artificial dissipation using the Laplacian would 
be A (Am) or in the semi-discrete case, 

v t = -Ii^P-'QaP^Qav. (15) 

In order to prove stability, equation (15) has to supplied with numerical 
boundary conditions. To investigate what numerical boundary conditions to 
use we examine the one-dimensional continuous counterpart of equation (15). 

u t = -h 4 u xxxx , 0 < x < 1 (16) 

Apply the energy method to equation (16). 

/ uu t dx = -h 4 uu xxxx dx = -h 4 uu xxx \l + h 4 u x u xx \\ - h 4 / u 2 xx dx (17) 

Jo Jo Jo 

In (17) it is seen that the boundary conditions h 4 u xxx = 0 and h 4 u x = 0 will 
result in a well posed problem. 

Next turn to the discrete equation (15) and apply the energy method to 
prove stability, 

\\v\\ 2 t = v T Pv t + vJPv = —h 4 (v T (A + DS)P~\A + DS)v + 

v t (A + DS) T P~\A + DS) t v) = 

—h 4 (2v T AP~ 1 Av + 2 v t AP~ 1 DSv + 2 v T DSP~ 1 {A + DS)v). (18) 

If equation (18) is compared with (17) we note that the boundary terms 
correspond to each other. In the discrete case we would choose, 

h 4 DSP~\A + DS)v = 0 , (19) 

h 4 DSv = 0 . (20) 

Equation (19) is precisely a discretisation of h 4 u xxx = 0 and (20) is the 
discretisation of h 4 u x = 0. To impose these numerical boundary conditions 
in a stable manner we again use the SAT technique (see [12, 13, 14]) and add 
the following terms to the right hand side of (15), 

h 4 P~ l DSP~ l (A + DS)(v - 0) + h 4 P~ l AP~ l DS(v - 0), (21) 

and obtain the artificial dissipation, 

v t = -h 4 P~ l AP~ l Av. (22) 
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To determine the size of (22) we consider the penalty terms (21). The sizes 
in (12) apply here also, which implies that the first penalty term in (21) is 
0(h 3 ), the second is 0(h 2 ) and equation (15) is 0(h A ). Since the action 
of the penalty terms is restricted to the vicinity of the boundary the global 
order of accuracy of (22) is 0(h 3 ). 

In fact, it is immediately realised that the artificial dissipation (22) might 
be a good choice because it is negative definite in the P-norm. However, the 
above derivation gives a more thorough understanding of the action of (22) 
and is also required to understand what order of accuracy that is obtained. 

The choice of numerical boundary conditions can also be understood di- 
rectly from the discretisation. The action of the penalty terms is to cancel 
the DS part of the second derivative discretisation. Thus, when this arti- 
ficial dissipation is used the Laplacian algorithm is run twice on the entire 
domain. The first time (v^b = 0 in (4) and the second time v n in (3) and 
(4) represents (Av) n , i.e. we have ((Au)jv)& = 0. This precisely corresponds 
to the boundary conditions u x = 0 and u xxx = 0. 

Remark In [11] it was shown that on grids different from equilateral poly- 
gons the approximation P~ l Q^v = const • Av + 0(1). Although, P - 1 Qa 
does not approximate the Laplacian on general grids it is always dissipative 
since A is negative semidefinite and has the correct size. Moreover, the dis- 
sipation is always a blend of derivatives, including the Laplacian. Hence, it 
has the same effect as the Laplacian and the additional dissipation will not 
affect a constant. 

Remark The fourth-order JST operator (see [9]) is built from a first-derivative 
approximation D\, a third-derivative approximation P> 3 and a scaling A (a:,) 
to become (DiXD^u. With this construction it is not possible to derive an 
energy estimate proving stability. 

5 Direction Dependent Artificial Dissipation 

In the previous section the artificial dissipation did not take directions into 
account. Sometimes, it is desirable to scale the artificial dissipation differ- 
ently in different space directions. In order to include that possibility, we 
will modify the approximation of the Laplacian. 

Following the analysis in [11] we consider, 

Pv t = Av, (23) 
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or at a specific point r,. 


v,M< = E 


neNi 


l 7 \ ^ dSin \ ^ dSi n \ ^ 

-dSi n — y ^ ^ ^ ^ di n V n H“ 

n£Ni Tm n£Ni Tm neN { 


where N t is the set of indices of neighbours to a point r t , implying that, 

\ ' dSin 

a ™ ~ ~ 2.^ ~r~ ’ 

^ AT 1 ni 
neNi 


dSl 


T ni 


if n € Ni. 


(24) 

(25) 

(26) 


Define a %3 = 0 whenever j 0 TV*. Further, a i3 is the (i, j) component of A. 
As is shown in [11], A has the following properties. 


-a r 


^y ^ U'lj — 'y ^ dm — 

i^zj neNi 

Qij — dji") 

an < 0. 


(27) 

(28) 
(29) 


These properties are equivalent to A being symmetric and negative semidefi- 
nite, which are necessary and sufficient for stability. That is, we may modify 
A such that the properties (27)-(29) are not violated and still obtain a stable 
approximation. In the case of artificial dissipation it will also be accurate 
(have the correct size) due to the scaling. 

Let e denote a unit vector in space in which direction the artificial dissi- 
pation should act and h tn the unit vector pointing at the direction r t — r n . 
We propose the following modification. 

Vi(v t )i = ^ — — — ds in \h in • e|. (30) 

n£Ni Tni 


In this case, 


(loin | ^ 

din = | m &\i 

Tni 

such that property (28) is fulfilled by definition. Also, it is easily seen that 
property (27) is satisfied (c.f [11]). Then (29) will also be fulfilled since cq n > 
0. The properties (27)- (29) are all satisfied, implying that this modification 
of the Laplacian is stable. 
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Equation (30) is just one example of a Laplacian based artificial dissipa- 
tion. A generalisation would be, 

Vi{v t )i = ^2 — — —ds in c ni . (31) 

_ at T ni 

neNi 

Stability will not be destroyed if c** = c. t] > 0 (c.f the above derivation 
and [11]). Thus, for example, it is possible to change the direction that the 
dissipation is acting at different parts of the grid without formally destroying 
stability. 

Remark With this observation in mind we would be able to choose c l3 = 
\ftij • e x + cf, where c is some constant, in the Cartesian case such that a 
nonzero dissipation is always obtained in the direction normal to x. 

To obtain a fourth-order dissipation we would apply (31) twice and ac- 
cording to (22) multiply it by h 4 where h is some grid size. The remaining 
issue is the choice of h. With a general unstructured grid the cell sizes may 
vary considerably in the domain and a better scaling should be used. How- 
ever, ri n is the size of an edge locally and we may use that when deriving the 
dissipation. 


( V t)i 


h 2 


-E 


V n Vi 


dSi n C n i ^ 


Tl^zNi 


- Yr 2 - 

V ^ n% 

1 neNi 


V n 


dsi n c n i 


1 

V 


Ek 

neNi 


Vi)T n idSi n C n 


(32) 


Again, conditions (27) -(29) are fulfilled since rj n = r n i > 0. 

The equations (22) and (32) constitute the final artificial dissipation. 
(Both are similar to the artificial dissipations used in [2, 15, 16]). That 
is, at each time step (32) is used once and yields an approximation of the 
Laplacian with homogeneous boundary conditions, (Am)*, at each grid point 
i in the domain. Equation (32) is then applied again on the grid function 
(Am)* instead of m* such that the fourth order dissipation is obtained. Note 
that with this approach, the edge based data structure is not destroyed, 
maintaining an efficient implementation of the finite volume scheme. 

Let us analyse the action of the proposed artificial dissipation, equation 
(32), on an equidistant Cartesian grid where the desired direction is e x . In 
this case we can order the solution in a matrix {%} where the first index 
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denotes space variation in the x-direction and the second index in the y- 
direction. Let h denote the grid spacing. In this case r t] = dsij = h and 
Vi = h 2 . We have for a point p, 

{Vt)pp = j^2 (('Up+i,p ^ p,p T Vp—i,p Xp,p)h |e® • e x \ + 

{vp,p+i v p : p “I - Vp^p—i v PjP )h \e y • c x | ) , 


or, 

ivt)pp = ^p+i,p + Vp—i^p, 

i.e. the standard second order undivided finite difference approximation of 
the second derivative in the x-direction. The Laplacian used twice will of 
course result in the standard fourth-order undivided finite difference approx- 
imation. 


5.1 Concluding Remarks on the Fourth Order Dissi- 
pation 

Consider, 

v t = -rfp-'AAP-'Av, (33) 

where A is a diagonal positive definite matrix. The energy method leads to, 

\\v\\ 2 T = -2v T AAP~ 1 Av = -h 4 2v T AA 1/2 P~ 1 A 1/2 Av, (34) 

where the right-hand side is a negative-definite quadratic form. Introduce 
A = AA 1 ! 2 . Equation (33) becomes, 

v t = - h i p- l A T P~ 1 Av , (35) 

A = AA 1 ! 2 can be interpreted as c VJ ^ c Tl in (31). To be more specific, 
c ni = c n for ah n E Ni and c lt > 0 can be chosen arbitrarily. This form of 
dissipation does not take directions into account but scales the dissipation 
only with respect to position in space. 


6 Numerical Computations 


6.1 Linear examples 


We will consider the two-dimensional advection equation. 


u t + au x + bu y 

Lu 

u(x,y,0) 


0, (—1 < x < 0, 0 < y < 1) = fl, 

g{x,y) (x,y)ed fi, 


(36) 
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where, 


1, (a, b) • n < 0 

0, (a, b) • n > 0 ’ 


(37) 


L = 


and n is the unit outward pointing normal. We will use the first derivative 
finite volume operators defined in [6] denoted P _ 1 Q X and P~ i Q y in the x- 
and ^-direction respectively. (Those are the standard scheme obtained with 
Green’s theorem.) In [6] equations such as (36), were proven stable with a 
weak implementation of the boundary conditions and we refer to that article 
for details. The second-order artificial dissipation operator is defined in (32), 
and repeated here for convenience, 


(D 2 v)i 


1 

V 


( Vn 

neNi 


^i^TnidSi 


(38) 


In the computations we use c n i = 1. To obtain a first-order dissipation we 
divide by r m , 


(.Diu)j — ^ ^ {pn Vi)dsi n . (39) 

* neNi 

Finally we define D±v = D 2 D 2 v. A semi-discretisation of (36) is, 

v t + aP _ 1 Q x v + bP~ l Q y v = e\Div + 64 D 4 U + BC, (40) 

where BC are penalty terms imposing the boundary conditions. Two test 
cases, computed on an unstructured triangular grid with 4357 nodes, are 
considered: 

1. ci = 0, e 4 = 0 or 1, random numbers as initial data, see Figure 2 . 

2. Ci = 0, e 4 = 0 or 5, sine function with random perturbation as initial 
data, see Figure 2. 

The results from Test 1 are displayed in Figure 3. Clearly, the dissipation 
damps the solution. Further, the non-dissipative computation converge very 
slowly to the steady state solution u = 0. Test case 2, Figure 4, also signi- 
fies the damping performed by the dissipation operator without altering the 
underlying solution. 
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(a) Random initial data. (b) Since function with random per- 

turbation. 

Figure 2: Two different initial data. 



(a) Without dissipation. 


(b) With fourth-order dissipation. 


Figure 3: Solution at t=l with random initial data. 
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(a) Without dissipation. 


(b) With fourth-order dissipation. 


Figure 4: Solution at t=0.5 with perturbed sine function as initial data. 


6.2 Nonlinear examples 


Consider Burgers’ equation, 



u{-l,y,t) 

u(x,y,0) 


0, (—1 < x < 0, 0 < y < 1) = e, 

u L , u(0,y,t) = u R , 

fi, 2 (x,y), 


where we use either, 


fi(x,y) 


ul — 1 , 

< -15a; - 12.5, 
k Ur = -0.5, 


-1 < a; < -0.9, 
-0.9 < x < -0.8, 
-0.8 < a; < -1, 


or 


h{x,y) 


Ul — 3, 

< -10a; -5, 
k Ur = 1, 


-1 < x < -0.8, 
-0.8 < x < -0.6, 
-0.6 < a; < -1, 


(41) 

(42) 


(43) 


(44) 


as initial data. We test 3 different cases where at most one of t\ and e 4 is 
nonzero: 

1. Initial data /i, unstructured fine mesh (11139 nodes). e\ = 1, e 4 = 0 or 
ei = 0, e 4 = 5. 



Figure 5: Solutions at t = 2 with different artificial dissipation. Test Case 1. 

2. Initial data /2, unstructured coarse mesh (2807 nodes) with either, 
ei = 0, 64 = 0 or, t\ = 1, 64 = 0 or ei = 0, 64 = 5. The shock will reach 
and pass through the boundary. 

3. Initial data / 2 , unstructured mesh with 11139 nodes with either, e\ = 
0, e 4 = 0 or, ei = 1, e 4 = 0 or e\ = 0, e 4 = 5. The shock will reach and 
pass through the boundary. 

In Fig 5 Test Case 1 is shown. We see that the first-order dissipation yields 
a smooth shock. The fourth-order dissipation is not able to damp all oscilla- 
tions at the shock but prevents them from spreading throughout the domain. 
In Fig 6, Test cases 2 and 3 are shown. On the coarse mesh all three options 
are stable. However, without dissipation the numerical solution has no ac- 
curacy and with dissipation the shock will move out through the boundary. 
The best performance is achieved with the first-order dissipation. On the 
fine mesh the non-dissipative even becomes unstable while the dissipative 
schemes manage to propagate the shock through the boundary. 

Remark We have tuned the amount of dissipation to be efficient in this case 
only. To be efficient in general flow calculations it would need a limiter that 
identifies shocks and yield a proper scaling of the dissipation. Nonetheless, 
our computations shows that the structure of the artificial dissipation scheme 
is efficient. 


7 Conclusions 

We have constructed and analysed first-, second- and fourth-order Lapla- 
cian based dissipation operators and proven them to be stable and accurate. 
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(a) Coarse mesh (b) Fine mesh 


Figure 6: Deviation from steady state solution as a function of t. Test Cases 
2 and 3. 

The first-order dissipation is suitable for shocks and the fourth-order dissipa- 
tion for non-physical oscillations. Computations corroborate the theoretical 
results. 
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